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Abstract 

This article deals with the spatio-temporal sensors deployment in order to maximize detection probability of an intelligent 
and randomly moving target in an area under surveillance. Our work is based on the rare events simulation framework. More 
precisely, we derive a novel stochastic optimization algorithm based on the generalized splitting method. This new approach 
offers promising results without any state-space discretization and can handle various types of constraints. 
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1 Introduction 

Consider a randomly moving object that tries to cross 
a closed area or to run away from a known position. We 
want to compute the best spatial and temporal deploy- 
ment of elementary search efforts in order to maximize 
the detection probability of this intelligent and randomly 
moving object. In our article, we call this object a target 
and the search efforts are assigned to sensors. Our work 
is closely linked to search theory [20], Stackelberg games 
[10] (which belong to game theory domain) and opera- 
tional research. It is also related to sensor network [7]. 

Stackbelberg games [16] could have been used to model 
our problem, however, due to limitations it was not pos- 
sible. In this two-person game, the first player, the leader^ 
commits to a strategy first. Then, the second player, the 
follower^ tries to optimize his reward, considering the ac- 
tion chosen by the leader. There are many security do- 
mains for which Stackbelberg games are appropriated. 
For example, the ARMOR system [15] has been used 
for years at the Los Angeles International Airport. AR- 
MOR casts a patrolling and monitoring problem as a 
Bayesian Stackelberg game (an extension of Stackelberg 
games). Moreover, theses games have potential applica- 
tions for network deployment and routing. Since we aim 
to deal with continuous temporal and spatial spaces, we 
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are faced with two difficulties. In our case, expressing the 
players' strategies as vector or a function would be dif- 
ficult. In addition, the number of strategies is expected 
to be quiet large or infinite. 

Search theory is another domain of game theory which 
deals with search efforts optimization. Most of the prob- 
lems in this field are used to be modeled in discrete space 
and time, and therefore, operational research methods 
are often chosen to solve these problems. Eagle [8] and 
Washburn [22] were the first to study the problem of the 
optimization of spatio-temporal search efforts. Using a 
branch and bound solver. Eagle proposed to maximize a 
bound on the detection probability rather than the prob- 
ability itself. Other approaches [1,3] only focus on spatial 
optimization without taking the target intelligence into 
account. Golen, meanwhile, recently [11] used a stochas- 
tic approach, more precisely, an evolutionary algorithm 
[2], but only optimize the sensors' discrete position. 

Our approach, based on the splitting framework [6], is 
different because it allows us to optimize both space and 
time and we do not have to discretize these spaces any- 
more. In the sequel, we will try to point out that the 
temporal aspect is the most difficult part of our prob- 
lem. This population-based optimization method can be 
linked to evolutionary algorithms. 

In this article, we first present and model our problem. 
Then, we give the expression of the detection proba- 
bility of interest and explain how we have applied the 
rare events simulation and the splitting algorithms to 
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our problem, resulting in the generalized splitting for re- 
search efforts scheduling (GSRES) algorithm. Lastly, we 
illustrate our method with the datum problem and give 
prospects. 

2 Problem presentation 

In order to maximize the detection probability, until time 
T, of a smart and reactive moving target whose trajec- 
tory is unknown (and depends in practice on random 
variables), we have to deploy both spatially and tem- 
porally a limited number of active sensors in an area Q 
(called the operational theater). Generally, 1] is a con- 
vex polygon with at least 3 vertices, but in most cases, 
it is a simple rectangle. As soon as a sensor is deployed, 
it is powered on and starts consuming energy ; it be- 
comes out of service when its battery is empty. In this 
article, we don't study the sensors' power consumption, 
so therefore we only assume that a battery is powerful 
enough to make the sensor emit a maximum number of 
times during a fixed period. Because sensors are active, 
we assume they can only detect a target when they send 
a signal. Since they are not autonomous, they only ping 
when a moving commanding station requests it. 

2.1 The target dynamic 

We are only given an a priori on the starting point of 
the target trajectory. This a priori is weak if the starting 
point is randomly sampled in the search space Q. On the 
contrary, a strong a priori means its initial position is 
sampled from a Gaussian pdf centered on the last seen 
position. 

We then define a feasible trajectory Yt G 3^ by a set of 
K points {oT K — 1 legs) composing it. So [17] 
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Fig. 1. Example of a leg- by- leg rectilinear trajectory (3 legs) 
Trajectory of a non-reactive and detected target 




Fig. 2. The target starts from yo. After a short time, it is 
detected by sensor 2 and does not change its course. The 
trajectory ends at yr- In red, the detection circle. In orange, 
the counter-detection circle. 



Yt = {yk}k=o...K-i where 
y-^ ||rfe+i -r/ell^ ^ Hr^-i - II2 ^ ^ (1) 

with Vk the target speed on the leg k. In this formula we 

note yk = [rx,k,Vx,k,ry^k,Vy,k]k=o...K-v where r^^k and 
ry^k define the target position and where Vx^k and Vy^k 
define the target velocity. 

Moreover, if the target detects a signal of a sensor but is 
to far from it, it is able to avoid it before being detected 
and it is also able to memorize its position. 

Basically, a target is instantly and surely detected if it 
enters the sensor's detection range. However, since the 
target is smart, if it comes close enough to a sensor which 
has just sent a signal, it detects it and learns all of its 
specifications. Thus, it may decide to avoid this threat or 



to come closer and start an avoidance later. In all cases, 
when the target is notified of the existence of a sensor, 
it changes its course before it enters the detection range. 
The target trajectory is then directly influenced by the 
partial knowledge of the target on the deployed search 
efforts fit{X). 

For instance, iJit{X) = means that until time t < 
the target has not detected any search effort yet. On the 
contrary, if lJit{X) = {1,2}, the target has learned about 
a sensor existence and is trying to avoid it (//^(X) = 1) 
or to run away from it after a detection {jit^X) = 2). 
This instant of detection by a sensor s is denoted by 

j-detect 

We also introduce the control vector u{P) which contains 
information about the change of course /3. At last, we 
can give an expression of the motion state equation of 
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Trajectory of an escaping target after a detection 




F{AtkiHh{X))) = 



Fig. 3. The target starts from yo. After a while, it is detected 
by an active sensor (sensor 1) and immediately escapes by 
following a radial course. 



Trajectory of avoiding target 




1 AtkiM^)) 
10 

1 AtkiM^)) 
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(3) 



Fig. 4. The target starts from yo. After a while, it detects 
an active sensor (sensor 2) and avoids it. 



More precisely, Atk is sarapled frora a truncated Gaus- 
sian law if /i/c(X) = {Atk ~ ^ft{l^At,cr%,/^At - 
ci-,IJ^At + ^ M). When a detection occurs {i.e. 

IJik{^) = {I52}), we end the current leg at tf^*^^* and 
add an extra one, so that K = i^T + 1, tk-\-2 = ^/c+i 
and t/e+i = t^^^^^^. Atk-\-i is then sampled from an uni- 
form law as Atk+i - l([tf ^^^S t/,+2]) (with tf^^^^ - 
U{[tk^tk-\-2]))' Whatever the value of /i/c(X), /3k is sam- 
pled from a truncated Gaussian law whose mean and 
bounds may vary depending on the last most threaten- 
ing sensor met. 

meanwhile, is sampled from a truncated Gaus- 
sian law {(3k ~ J\ft{lJi^,cF'^^,IJi^ - a.fi^ ^ a), a e M) 
when /i/c(X) = {0,1}, and from an uniform law when 
/i/c(X) = 2. More precisely, when /i/e(X) = 0, the mean 
of the Gaussian law is equal to the initial course of 
the trajectory or to its last course Pk-i^ depending 
on the chosen dynamic. When iJLk{^) = 1, is directly 
defined by the last most threatening sensor met. Lastly, 
if the target is detected, i.e. iJik{^) = 2, the new course 
of the target is defined by its previous course and the 
position of the detecting sensor. 

2.2 The solution 



Remember that we have a set of P sensors Si that we 
may spatially and temporally deploy in our operational 
theatre Vtx [0, T] in order to detect a smart and reactive 
target. A sensor is only operational for a limited amount 
of time {e.g. its battery life) and can therefore send a lim- 
ited number of pings {Emax times maximum). Moreover, 
sensors are not autonomous and they are only able to 
send a signal if they are given the order by a command- 
ing station. This is denoted by the visibility parameter 
until time t < T, (^t{X). 

Before going further, we need to explain the following 
points: 



our target: 



Vk+i = F{Atk{fik{X)))yk + Ut{l3k{fik{X))) . (2) 



Here, Atk and Pk are two random variables and 
F{Atk{/J^k{X))) is a state transition matrix associated 
with a rectilinear motion: 



• all the sensors must be in the search space 

• sensors must have been set up before they can be ac- 
tivated, 

• instants of activation must be in [0; T]. 

Denote a solution by X G A' (the set of feasible solu- 
tions), so we can define: 

X = {Si(Ti)}i=i,...,p with Ti = {ti^i, . . . , tij}j = i^,,,^npi- 

(4) 
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Si corresponds to the sensor i {i G {1, . . . , P}) position, 
while Ti is a vector containing its npi instants of activa- 
tion (in [0, T]). We also denote by C all the spatial and 
temporal constraints on the feasible solutions. 



3 Evaluating the detection probability 

We want to maximize the detection probability of a tar- 
get until time T. This quantity is denoted by St{X) and 
is given by the following equation: 



St{X) ^ f f {YT\^T{X);C)p{YT\fiT{X);C) dYr 



(5) 



Here, f {Yt\^t{X)]C) is a cookie-cutter cost function 
that takes the value 1 if the studied trajectory Yt sat- 
isfies some defined criteria (such as a number of detec- 
tions and a number of avoidances) and otherwise, i.e. 
f{YT\MXy,C) = l{v,GA(x,c)} where A{X,C) is the 
set of trajectories which are detected by the solution X 
and which respect the constraints C. It depends on the 
visibility of the solution ipT{X). p (YT\fiT{X);C) is the 
conditional pdf used to generate the target trajectories 
and depends on the target intelligence iit{X). As this 
is not the goal of this article, we do not give any more 
information on how we have implemented this cost func- 
tion. Unfortunately, St{X) is an integral with respect 
to the probability distribution of the (random, solution- 
dependant) target trajectory Y and its analytical ex- 
pression is not available. A first approach should be to 
use the drude Monte Carlo method to obtain an unbi- 
ased estimator of St{X)^ St{X): 



N 



St{X) = i^/(Y^|^r(X);C), where 



1=1 



(6) 



YS ^ p{YtW{X)-C). 



The trajectories 1^ are recursively generated using the 
motion state equation (2). To be concrete, we generate 
a large number of feasible trajectories 1^, i = 1, . . . , TV 
and evaluate / (l^ | ifT (^) ; C) . 

Note that the relative error associated with St{X) given 
by the CMC estimator is 



REcMc{Sr{X)) = -^== 



(7) 



Also remark that the smaller the probability to estimate 
is, the larger the relative error is. To reduce this error, we 
have to increase the number of trajectories N . Knowing 



the probability we are meeting are above 10 ^ (if they 
were below, planning would be useless), we have chosen 
TV > 50000. 

The optimal solution we search, denoted by X* is de- 
fined by this formula: 

^ argmax^^;,{5T(X)}, XeX. (8) 



We also define 7^ = St{X'^) as the maximum detection 
probability we can reach. Let X"^ = arg max{S'T(X)} be 
the solution obtained when maximizing the approached 
probability 7^ = 5't(X'''). For the sequel we focus on 
computing 7"^ and X"^. Indeed, if TV is large enough, we 
can expect that these are good approximations of 7"*" 
and X^ respectively. A suflScient condition for this to 
hold is that St{X) converges to St{X) when N ^ 00 
uniformly in X G A'. 

Splitting theory is based on the observation that maxi- 
mizing St{X) (in practice St{X)) is similar 

• to estimating the probability 

£(7) = nMx) > 7) = / 

(9) 

where g(X, C) is some positive probability density on 
A', with the idea that this probability decreases and 
reaches zero when 7 increases toward the (unknown) 
value 7^ and is identically zero beyond this value, 
hence the characterization 



7t = inf {7 > : ^(7) = 0} , 
and simultaneously to sampling the set 



(10) 



X^ = {X^X'. St{X) > 7} C A' , (11) 

with the idea that this set decreases toward X''' when 7 
increases toward the (unknown) value 7^ and reduces 
to beyond this value. 

However, when 7 goes to 7''', the event {St{X) > 7} 
becomes rarer and rarer, and consequently the CMC 
estimator 

1 ^ 

^(7) = ^ E ^{S.(X.)>7} ^^^^^ - ^) (12) 

i=l 

has a relative error 



(13) 
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that increases to infinity. In order to reduce this relative 
error, we should increase the sample size C, but then, we 
would be faced with a computational explosion. More- 
over, when 7 goes to 7''', it becomes harder and harder 
to produce samples from C) that would be close to 



To address this problem, a technique called generalized 
splitting [5] and derived from Diaconis, Holmes and Ross 
researches [9] on MCMC (Markov chain Monte Carlo) 
allows us to compute ^(7) in an easier and more precise 
way. For an optimization problem, we will find at least 
one solution that maximizes our criteria among all the 
solutions sampled to compute our probability of interest. 



If we define a sequence of increasing thresholds 7, such 
that 7o > 7i > • • • > 7l (with ^ 7"^), we can rewrite 
^(7) as the following product of conditional probabili- 
ties: 



importance sampling function, ^(7) can be rewritten: 



1^ 

£(7) = P, (St{X) > 70) {St{X) > 11) 



L 
1=1 



(16) 



With a judicious choice or a fair estimation of the {7/} 
sequence, the event {St{X) > 7/} is no longer a rare 
event (generally, q G [10~^, 10~^]) under the distribu- 
tion gi_-^{X ^ji-i^C) and therefore the q quantities can 
now be well approximated through a CMC estimator. 
Hence, a CMC estimator of ^(7) is: 



(17) 



L 

n 

1=1 



^(7) = Fg (St{X) > 70) (St{X) > ^i\St{X) > 7i_i 



where q = ^ Ylf=i ^ 
^f_i(X;7z_i,C). 



{5t(X,)>7J 



and where X,- 



CO 



where 



(14) 



ci = P, {St{X) > ^i\St{X) > 7^_i) 
P, {St{X)>7i.St{X)>7i-i) 



-L 
-L 



dX 



(15) 

and where the importance sampling density [18] 
gU{X,^i-„C) = is precisely 

the conditional density of X, given that St{X) > 7/_i. 
Remark that the support of this density gi_-^{X] 7/_i , C) 

is precisely the set {X G X : St{X) > 7/-1}. 



If we know how to draw independent and identically dis- 
tributed random variables Xi over G X from this 



4 The splitting algorithm 

We will now describe the general splitting algorithm ap- 
plied to our optimization problem. 

Generate a population = {Xi, . . . ^ q{X]C) 

of C feasible solutions (that respect all the constraints 
in C) and initialize an iteration counter / = 1. Evaluate 
the scores of the solutions Sq = {SriXi)} and sort So 
in decreasing order such that S't(^j(i)) > 5'T(^j(2)) ^ 
... > 5't(Xj(c)). We obtain 70 = 5't(Xj(Co)) with Co = 
[pC\. Set = {^}^=l,...,Co C such that SriXi) > 
7o. Remark that Xi ~ ^0 (-^5 lo^C) for i = 1, . . . , Cq- 

At iteration / of the algorithm, contains only p% 

samples of the initial population whose score is 

above the current threshold 7/-1. Since we want to 
work with a constant number C of solutions, we have 
to complete this set. This is done in two steps. First, 
we use the bootstrap method or the ADAM cloning 
method to repopulate our set of solutions. The boot- 
strap method consists of sampling uniformly with re- 
placement C times from the population ADAM 
cloning method, meanwhile, consist of making -^-^ + 

Bi{i = 1, . . . , copies of each sample. Here each 

Bi, . . . , Bi-i are Ber {1/2) random variables conditional 



on J2 



Ci- 



B, 



C mod Ci-i and . . . , Bci_^] is a 



binary vector with joint pdf P(5i = 61, ... , Bci_ 

^Ci-^) = -^%^f^l{6i+...+6c,_i=r}, bi G {0,1}, where 
r = C mod C/_i. Unfortunately, the samples of the 
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completed set denoted by x\^^^^^^^ are identically dis- 
tributed but not independent. 

To address this problem, we apply a random Gibbs sam- 
pier 7r,_i(X|X,_i;C) = ^ «,_i(X|Xi;C) to 

each sample Xi of x^^^^^^^^^ to obtain Xi = {Xi} such 
that Xi ^ gf_-^^{X;ji_i,C) for i = 1, . . . , C Here, 
is the transition density of a Markov chain starting from 
^hoot/cion with stationary pdf g^_i. For our prob- 
lem, we have the transition density defined by: 



(X\X^ = E fi , (18) 

j=l r=l 

where XI denotes the component r of a solution and 

all the components of X^ but the r one. The Xj are 
the probabilities of updating one component at a time, 
given that Xj = 1 and the ruj are the conditional pdf 
associated to the 6 moves (described below). We could 
also use a systematic Gibbs sampler, whose transition 
density is defined by 



1^1 



(19) 



r=l 



Notice that we have chosen a random Gibbs sampler with 
bi random updates of components of a solution Xi , while 
with a systematic Gibbs sampler, we would have updated 
all the components of Xi in a fixed order. The number 
of random updates bi varies during the simulation in 
this way: bi = bo -\- al where a G W^- For the first 
iterations, bi < P and therefore this approach is faster 
than a systematic Gibbs sampler. On the contrary, when 
/ is close to L, bi > P. Thus, we do more updates than 
a systematic Gibbs sampler would do but we maintain 
more diversity in our solutions. 

Since we do not know how to update a solution in a 
way it still satisfies the constraints C, we first recur- 
sively propagate the modifications starting from the sen- 
sor/activation we have modified in the sequence of acti- 
vations. Then we check its feasibility, that is, if it respects 
all the spatial and temporal constraints C. We apply ac- 
ceptance-rejection (for a limited times) to each updated 
component until we find a feasible solution. Considering 
that the cost function St also verifies the consistency 



of a solution, an updated solution Xi from A!- 
then accepted with probability ^^St(x-)>ji i} 
forth, all the Xi should be iid. 



boot /cl on 



. Hence- 



In the same way as before, evaluate the scores Si = 
{SriXi)}, Xi G P^i and sort this set in decreasing order 
such that Sri^jii)) > St{Xj^2)) > •.• > M^jiC))- 



We obtain % = ST{Xj(^Ci)) with Ci = [pC\. Deduct 
that Xi = ^j(i) = diTgmdix{ST{Xi)}, i = 

and 7/ = St{Xi) respectively the best solution at iter- 
ation and its score at iteration 1. In addition, denote by 
70:^ = niax{7/,7o:/_i} and Xq.j = Xi if 7^ > 7o:/_i or 
Xo:Z = Xo:/-i otherwise, respectively the best detection 
probability and its associated solution ever met until it- 
eration /. 

If one of the stopping criteria is satisfied, stop the al- 
gorithm and give Xq.j as an estimator of the optimal 
solution. 

In practice, the threshold 7/ is the p quantile of the sam- 
ple Xi. 

4-1 The Gibbs sampler moves 

Before we go further, let us introduce and recall a few 
notations. Si = [si^; Si^] denotes the ith sensor position 
(and more generally the ith sensor) , P is the number of 
sensors in the current solution, Pmax is the maximum 
number of sensors, npi stands for the activations' number 
for sensor i while and respectively are the 

instants of activation of sensor i and the set of activation 
times associated with sensor i. Also denote by t^. the 
set up duration of the sensor Si. Remark that a sensor 
whose instants of activation are negative is considered as 
disabled. Consequently, deleting an instant of activation 
consists of assigning a negative value to this instant. 
Removing a sensor is then equivalent to deleting all of 
its instants of activation and ignoring it. Now we give 
more details on the 6 moves of our Gibbs sampler and 
the associated conditional pdf. 

(1) Add a sensor. The conditional pdf mi can be de- 
fined in two steps. Sample a position Sp_^^ from 
U{fl',C) for the new sensor. Then draw its first in- 
stant of activation tp^^^ 1 ^ ^{[^sp+i^T]). mi is 
proportional to (up to a normalization constraint) : 



mi 



(20) 



(2) Add one instant of activation. First, choose 
a sensor randomly i.e. draw j uniformly in 
{!,..., P}. If npj < npmax then draw ^ 
ZY([tj,i, T]). 7712 is proportional to: 



7712 



U{{1,...,P}) W([i,,i,T]) l{s,(x,)>^,_,} . 



(21) 



(3) Remove a sensor. To apply the move 7713, we ap- 
ply the following steps : choose a sensor randomly. 
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i.e. draw j uniformly in {1, ... , P}. Then delete all 
of its instants of activation and mark it as disabled. 
So we have 



7Z-l} 



(22) 

(4) Remove an instant of activation. Choose a sen- 
sor randomly i.e. draw j uniformly in {1, . . . , P}. 
We assume that npj > 1. Choose an instant of ac- 
tivation tj^k^ draw k uniformly in {2, ... , npj}. 
Delete t'-j^. The conditional pdf 777.4 is defined as 
below: 



m4 



{Xi\xri) cx 
U{{1, . . . , P}) W({1, . . . , np,}) l{s^(x,)>^,_o 



(23) 

(5) Move a sensor. Select a sensor Sj randomly, i.e. 
draw j uniformly in {1, . . . , P}. Then, draw s'^ ^ 

J2k=i^k JV{sj,T,l) with Yfk=i^k = 1. Notice 
that the weights Wk may evolve during the opti- 
mization in order to promote one of the move ver- 
sus the other. For this mixture of two Gaussian pdf 
the covariance of the first Gaussian defines a small 
move while the covariance of the second Gaussian 
defines a larger move. Here is the formula of this 
conditional pdf: 



W({1,...,P}) Efc=i«'fcA/'(«^,S^)l{5T(x,)>7,-i} • 

(24) 

(6) Swap two sensors. If we assume there are at least 
2 active sensors, select two sensors Sk and Sr with 
k uniformly drawn in {1, . . . , P} and r uniformly 
drawn in {1, . . . , P} \ {k}. For all /c = 2, . . . , npk^ 
delete j and f^j for all r = 2, ... , npr- Swap their 
first instant of activation: ^ and ^ . The condi- 
tional bivariate pdf ttiq is then defined by: 

me (x^XriX-^'-^) oc 

Z^({1, . . . , P}) Z^({1, . . . , P} \ l{5,(x,)>^,_,} , 

(25) 

where e {(Xf , )}. 

^.^ T/ie GSRES algorithm 

Algorithm 1 (GSRES) Given parameter p, sample 
number C and number of hum-in iterations hi of the Gihhs 
sampler, follow the forthcoming steps: 

1, Initialization. Set a counter 1 = 1. Generate C 
feasible solutions {X^}, i = 1, . . . , C and denote Aq the 
set containing them. Note that X^ ^ g(X;C). Evaluate 
scores Sq = {SriXi)} and sort in decreasing order Sq 



such that ST{Xj(^i)) > 5't(X^(2)) > ••• > 5't(X^(c))- 
We obtain 70 = S'T(-X^j(Co)) '^^^^ ^0 = [pC\. Define 

Xq = Xo:0 = ^j(l); 7/ = 70:0 = 5'T(^j(l)). 

2, Selection. Let ^i-i = {Xi, ...,Xc^_^} 6e ^/le subset 
of the population {X I, Xc} for which St (Xi) > 7^-1. 

contains p% of the population. Notice that Xi ^ 
9i-i{X;^i_i,C) fori = 

3. Repopulation. Apply one of these methods: 

• Bootstraping: sample uniformly with replacement C 
times from the population ^i-i to define the temporary 
set of C solutions 



^hoot 



ADAM Cloning: make -]- Bi{i = 1,...,C/) 

copies of each population sample Xi-i. Here each 
Pi,...,Pq are Per(l/2) random variables condi- 
tional on Yl^Li Bi = C mod C/. We then define the 
temporary set of C solutions Af/ 



don 



4' Gihhs sampler. Apply a random Gibbs sampler 
^,_i(X|X,_i;C) = ^^f-^„,_i(X|X,;C) with h 
burn-in iterations and the transition density ki-i to each 
sample qJ^ p^^^^^/^^^^ f^g^^ section to obtain Xi = 
{Xi} such that Xi ^ gf_^{X;%-i,C) for i = 1, . . . ,C . 
Notice that the Xi^i = 1, . . . , C should be approximately 
independent and identically distributed. 

5. Estimation. Evaluate scores Si = {ST{Xi)}, Xi e 
Xi. Sort in decreasing order Si such that ^ 
ST{Xj(^2)) > ••• > SriXji^c))' obtain 7/ = 
SriXji^Ci)) '^'^th Ci = [pC\. Deduct that X/ = -X^j(i); 
7/ = 5t(X^(i)); Xo:i = Xi ifji > 7o:/-i; else Xq-.i = 
Xo:^-i and%:i = max{7/, 7o./_i}. 

6, Stopping condition. If one of the stopping condi- 
tion is reached, stop the algorithm and give Xq.j as an 
estimator of the optimal solution. Else I = I -\- 1 and go 
back to step 2. 

For our problem, we implement a customized version 
of the splitting method. To begin the computation, we 
generate an initial pool of feasible solutions with q{.;C). 
Since a solution and the carrier trajectory are closely 
linked, we use our trajectory generator to obtain a pool 
of initial solutions that respect the whole constraints set 
C. 

Lastly, to ensure our algorithm will not converge and 
stay into a local extremum, we developed a simple heuris- 
tic. If the current maximum score and the current thresh- 
old do not increase for a chosen number of times, we au- 
tomatically reduce the value of the threshold. Through 
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the decrease of the threshold, we start again to accept the 
feasible solutions generated by the moves and therefore, 
we reintroduce some diversity in the pool of solutions. 

5 Illustrative example 

The first result we present here concerns a scenario in 
which a target is running away from the position it has 
just been detected. Its initial position is drawn from 
a Gaussian law centered on 1^/2 and with a variance 
^target - Morcovcr, the target is supposed to be smart and 
reactive and therefore, while it is running away, it tries 
to avoid being detected another time. Considering that 
the research starts with a delay of t^^^ which represent 
the time of arrival of the hunter, we aim to maximize the 
chances to detect the target during the time T. We use 
Pmax = 10 sensors that are able to ping only once. For 
this simulation, we use C = 800 solutions, N = 70000 
trajectories, bo = 2, bi = bo -\- 0.2 I and decide to keep 
10% of elites {p = 0.1). We also let the algorithm per- 
form up to 50 iterations. 

Because our algorithm is not able yet to adjust the num- 
ber of sensors considering the cost of their deployment, 
we have chosen to work with a constant number of sen- 
sors. However, we have allowed the removal of a sensor if 
it is directly followed by an addition of a new sensor. We 
have used two of the six moves we have defined above : 
move a sensor and a combination of removing a sensor 
followed by the addition of a new sensor. The probabil- 
ity of both moves to occur is 0.5. 

In the best solution we obtain, the sensors position and 
activation describe a spiral. This result, illustrated in 
the figure 5, is related to the studies of Washburn [21] 
and Son [19] for an only-spatial optimization case, i.e., 
when the target is not able to avoid the sensor ("myope" 
case). In this context, the best spatial sensor deployment 
designs an Archimedean spiral. 

In the sequel, we plot the optimization evolution be- 
haviour versus iterations joa = StO^o-.i) and E[5't(X^)] 
which represents the mean score of the current popula- 
tion (figure 6). Both are smoothly increasing in a loga- 
rithmic way. On the figure 7 we see that the support of 
scores pdf, much large at the beginning (/ = 0), becomes 
thinner and thinner and converges toward a Dirac pdf as 
the optimization is conducted. Moreover, joa increases 
and the standard deviation of the distribution decreases. 
Once the optimization is over, we observe that the detec- 
tion probability reaches 0.9406 whereas when / = 0, the 

best score is below 0.25 and the mean scores E St{X) 

is equal to 0.0187 with a large standard deviation. This 
gap illustrates the efficiency of our approach. 

Figure 8, visualizes the result of two other simulations 
for wich the planification is not as optimal as for the pre- 



Sensors position and order of their first activation 




Fig. 5. Graphic of X' : position and activation order of the 
10 sensors. 

,Sr(X():/) and E[5T(Xi)] evolution with C = 800, N = 70000, L = 50 




Fig. 6. In blue StC^o-.i) and in red E[5't(X0] with C = 800, 
N = 70000, L = 50. 

Solutions distribution for 4 iterations of algorithm with C = 800, N = 70000, L = 50 



Sr(X,), {XJi.!,... 

§t(xo, {X0..1,... 
ST(Xi), {X0..1,... 



Fig. 7. Scores densities support versus iterations GSRES. In 
blue / = 0, in green / = 5, in red / = 10 and in black / = 50. 

vious example. Nevertheless, the detection probabilities 
associated with these deployments are above 0.93. 

To compare our best solution (see figure 5) with Son's 
solution, we have computed the best track spacing of 
the Archimedean spiral (denoted by TS) for our sce- 
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Fig. 8. Graphic of \ position and activation order of the 
10 sensors for two other sub-optimal solutions. 

nario (with CMC method for detection probability eval- 
uation), considering the target has a myope behaviour. 
According to Son's researches, the best track spacing is 
computed thanks to the following equation: 

TS = msix{2R,mm{aTSray^TS'' ^ P}} . (26) 

TSray IS Considered as a good track spacing for a ran- 
dom tour target ^ that is approximately normally dis- 
tributed, S'^ is the largest track spacing (so-called the 
"furthest-on-disk" solution) and R is the sensor detec- 
tion radius. In [19], a and j3 have been obtained by us- 
ing a quasi Monte Carlo framework, more precisely, a 
nearly orthogonal latin hypercube (NOLH) method. Af- 
ter that, we have fitted an Archimedean spiral and a log- 
arithmic spiral to our solution (see [14]). Figure 9 shows 
the comparison between the myope Archimedean spiral, 
our solution, the fitted Archimedean spiral and the fitted 
logarithmic spiral. We remark the sensors position and 



^ Assuming that the course change frequency is large enough 
versus ^ and that the carrier speed is much larger than the 
target speed. 



Datum optimal solution and fitted Archimedean and logarithmic spirals 

Optimal solution 

Fitted Archimedean spiral 

Fitted logarithmic spiral 

+ Experimental solution 




Fig. 9. Datum optimal solution and fitted Archimedean spiral 



Detections per sensor 




Fig. 10. Sensors detection rate for the best solution X' . 

activation rather describe a logarithmic spiral than an 
Archimedean spiral. This may be explained by the dif- 
ferences between the target dynamic models and by the 
fact that we don't use a single and always active sensor, 
but 10 sensors instead. 

Figure 10 shows that the first sensor deployed detects 
most of the targets (TV = 70000). 

The following pictures describe how our algorithm 
works. Figure 11 shows the spatial distribution of the 
C = 800 solutions at initialization. We remark that for 
the first four sensors, the spatial distribution densities 
follow a Rayleigh distribution. As the carrier trajectory 
bounces in the search space for the next sensors, the 
spatial distribution densities converge toward a uniform 
distribution. 

In figure 12 we notice that the first sensor is almost 
positioned. For the other sensors, some saddle points 
emerge. 
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Fig. 13. Spatial sensors' pdf at final iteration 

At final iteration (see figure 13), the sensors are posi- 
tioned and all the spatial distribution follow a unimodal 
distribution. 

During the Gibbs sampler's step of the optimization, we 
reject a new solution if it is not feasible or if its score 
is below 7/_i at iteration /. Otherwise, we consider that 
the move is accepted. The two following graphics (fig- 
ures 17 and 18) show the evolution of the two moves' ac- 
ceptance during the optimization. The acceptance rates 
decrease while the optimization is conducted. The two 
peaks correspond to the threshold's decrease. 

The next ones (figures 19 and 20) focus on the move's re- 
jections. More precisely, the blue curves describe the rate 
of non-feasible solutions generated by the move, and the 



Fig. 16. Temporal sensors' pdf of sensors at final iteration 

green curves describe the rate of solutions which score 
is below the threshold ji-i conditional on the fact they 
are feasible. Here again, the sharp drops of the rejections 
designed by the green curves correspond to the thresh- 
old decrease. 

6 Conclusion and prospects 

In this paper, we presented a method to compute the 
best deployment of sensors in order to optimize the de- 
tection probability of a smart and reactive target. Unlike 
existing works that use discrete constraints or only fo- 
cus on the optimization of one aspect at a time, we pur- 
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Fig. 17. Acceptance rate for move "delete a sensor" 




Fig. 18. Acceptance rate for move "move a sensor" 




Fig. 19. Rejection rates for the move "delete a sensor" 

posed to optimize both space and time with continuous 
constraints. This was made possible through the use of 
a novel stochastic optimization algorithm based on the 
generalized splitting method. 

We tested our algorithm with the datum search prob- 
lem and the results we obtained were very satisfying. In 
addition, these results were confirmed by existing works 
on much simplier cases. This demonstrate the efficiency 
of our dedicated Gibbs sampler and the splitting frame- 
work for this type of problem. 

To reduce the rejection rate and consequently reduce the 
computing time, we will try to improve the efficiency of 



Fig. 20. Rejection rates for the move "move a sensor" 

our moves. In the same time, we should focus on coop- 
eration between sensors {i.e. multistatic case) and use a 
more complex model for sensors, signal propagation and 
detections {i.e. not cookie-cutter sensors). We also aim 
to take the cost of each solution into account and com- 
pute a real multiobjective optimization. To achieve this, 
we are currently investigating two different ways. The 
first one would be based on a Pareto-ranking algorithm 
[4] and the second would use Choquet integral for the 
aggregation of the user preferences [12,13]. 
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